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ABSTRACT 


Improving  the  current  Middle  East  wave  speed  with  full  waveforms  required  confidence  in  sources  and  recordings, 
along  with  a  methodology  to  iteratively  improve  the  model  and  reducing  its  minimum  period.  Recordings  of  seismic 
waves  traversing  the  region  from  Tibet  to  the  Red  Sea  and  their  mismatch  compare  to  synthetics  from  the  current 
iteration  model  are  the  principal  metric  in  improving  the  current  wave  speed  model.  To  avoid  a  mapping  of  source 
errors  into  the  new  wave  speed  model,  a  rigorous  characterization  of  each  source  within  the  current  wave  speed 
model  was  undertaken.  Source  depths  and  paths  near  nodal  planes  are  error  prone  as  small  changes  may  affect  the 
resulting  wave  field.  After  the  sources  were  evaluated,  regions  requiring  refinement  were  highlighted  using  adjoint 
tomography  methods  based  on  spectral  element  simulations  (Komatitsch  and  Tromp,  1999). 

We  reinterpreted  a  large,  well  recorded  subset  of  201  events  (1997-2007)  through  a  direct  comparison  between  data 
and  synthetics  based  upon  a  centroid  moment  tensor  inversion  (Liu  et  al.,  (2004).  Initial  evaluations  were  done  using 
a  ID  reference  model  (Dziewonski  and  Anderson,  (1981)  at  periods  greater  than  80  seconds  and  a  more  stringent 
evaluation  was  done  for  three-dimensional  (3D)  model  (Kustowski  et  al.,  2008)  at  periods  of  25  seconds  and  longer. 
Final  source  reinterpretations  within  the  3D  model  define  a  source  database  and  the  initial  starting  point  for  the 
adjoint  tomography.  Transitioning  from  a  ID  to  3D  wave  speed  model  shows  dramatic  improvements  when 
comparisons  are  done  at  shorter  periods  (25  s)  and  even  at  longer  periods  (80  s).  Synthetics  from  the  ID  model  were 
created  through  mode  summations  while  those  from  the  3D  simulations  were  created  using  the  spectral  element 
method. 

Finally,  updates  of  the  wave  speed  model  were  accomplished  using  adjoint  tomography.  Initial  inversion  iterations 
have  guided  the  measurements,  model  smoothing  and  filtering  to  produce  the  accurate  and  relevant  improvements  to 
the  initial  model.  Initial  attempts  to  update  the  wave  speed  model  were  hampered  by  the  strong  anisotropy  in  the 
mantle  causing  an  unavoidable  mismatch  between  Rayleigh  and  Love  waves  when  using  an  isotropic  mantle 
parameterization.  Relaxing  constraints  on  the  mantle  wave  speeds  to  allow  for  transverse  isotropy  with  a  vertical 
symmetry  axis  allowed  the  fitting  of  both  surface  waves  and  thus  continued  improvement  of  the  wave  speed  model 
and  inclusion  of  even  shorter  periods. 

Event  kernels  are  computed  “quickly”  at  the  n-th  iteration  and  invert  for  a  new  (n+l)-th  wave  speed  model  of  the 
Middle  East.  As  demonstrated  with  previous  adjoint  tomography  experiments  (Tape  et  al.,  (2009),  each  iteration 
improves  the  model  progressively,  and  relies  on  successive  iterations  to  reduce  variance  and  improve  the  fitting  of 
data  to  synthetics.  Exploiting  the  source  database  with  multiple  adjoint  inversions  at  shorter  and  shorter  periods  will 
refine  the  wave  speed  structures  of  the  Middle  East.  Inversion  results  demonstrate  that  the  iterative  nature  of  the 
adjoint  tomography  improves  the  travel  time  variations  between  synthetics  and  data  primarily,  with  less  of  an 
improvement  to  the  waveform  amplitudes.  Iterative  improvements  also  significantly  increase  anomaly  strength 
while  sharpening  the  anomaly  edges  to  create  stronger  and  more  pronounced  tectonic  structures.  The  results 
presented  here,  while  accurate  at  intermediate  periods,  require  the  addition  of  attenuation  tomography  and  a 
transverse  isotropy  without  a  vertical  symmetry  axis  to  further  reduce  the  minimum  period  towards  travel  time 
tomography  models. 
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OBJECTIVES 


Improved  3D  wave  speeds  of  the  greater  Middle  East,  from  the  Turkish  Plateau  to  the  eastern  edge  of  Tibet 
(Figure  1),  enhances  our  ability  to  discriminate  between  natural  and  man-made  events,  locate  these  events,  identify 
source  depths,  and  determine  magnitudes  using  waveform  based  methodologies.  Current  wave  speed  models  of  the 
Middle  East  are  improved  through  an  adjoint  tomography  method  (Tromp  et  al.,  2005)  by  a  comparison  of  seismic 
phase  and  amplitude  data.  An  initial  step  towards  the  adjoint  tomography  model  requires  the  development  of  a 
database  of  relocated  and  well-characterized  set  of  sources  and  waveforms.  Events  were  re-inverted  in  ID  and  3D 
wave  speed  models  and  agree  fairly  well  with  currently  available  solutions.  Use  of  a  3D  model  improves  the 
waveform  fit  between  data  and  synthetics  and  reduces  the  overall  error  of  the  CMT  solution.  The  adjoint 
tomography  method  (Tromp  et  al.,  2005)  uses  full  seismic  waveforms  as  a  measure  of  misfit  of  the  current  model 
iteration.  Differences  between  data  and  synthetics  are  used  to  create  adjoint  sources  and  generate  sensitivity  kernels 
required  to  update  the  wave  speed  model.  The  final  round  of  iterations  amplifies  and  sharpens  wave  speed  anomalies 
in  the  Middle  East  while  increasing  the  predictive  capabilities  at  periods  of  15  seconds  and  longer.  The  adjoint 
tomography  model  was  built  on  the  foundation  of  the  seismic  waveform  database  and  finite-frequency  wave 
propagation,  and  will  be  distributed  to  the  community. 

RESEARCH  ACCOMPLISHED 


Event  Characterization 

Approximately  200  events  in  the  Middle  East  were  re-characterized  using  the  CMT  inversion  methodology  of  Liu  et 
al.  (2004)  and  the  ID  Preliminary  Reference  Earth  Model  (PREM)  (Dziewonski  and  Anderson,  1981)  at  periods  of 
25  /  60  seconds  and  longer  and  within  a  3D  wave  speed  model  (S2.9EA)  (Kustowski  et  al.,  2008).  A  map  of  all 
available  stations  and  events  are  shown  in  Figure  1.  Initially,  the  procedure  of  processing  large  amounts  of  data, 
comparing  these  to  synthetics,  and  reevaluating  the  source  parameters  and  locations  needed  to  be  assessed  and 
streamlined.  This  procedure,  accomplished  through  the  use  of  a  ID  wave  speed  model  and  long  period  modeling, 
needed  to  be  straightforward  and  nimble  enough  to  avoid  problematic  areas,  as  the  initial  3D  event  re-evaluation 
requires  more  demanding  computational  efforts.  Problematic  events  and  stations  were  subsequently  screened  out 
within  this  procedure  before  building  the  final  waveform  data  set  was  implemented  in  the  adjoint  inversion. 
Problematic  stations  and  waveforms  with  dropouts  and  poorly  characterized  amplitude  responses  can  negatively 
influence  a  CMT  inversion,  and  these  stations  were  removed  before  any  3D  CMT  inversions  or  adjoint  inversions 
were  performed. 
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Figure  1.  Map  of  the  available  sources  and  stations  in 
the  Middle  East  for  the  adjoint  tomography 
inversion. 


Solutions  using  the  3D  wave  speed  model  at  periods 
of  25  seconds  and  longer  are  shown  in  Figure  2. 
Faulting  parameters  agree  well  with  those  from  the 
Global  Centroid  Moment  Tensor  (CMT)  Project. 
Solutions  match  reasonably  well  when  using  either 
ID  or  3D  wave  speed  models,  but  use  of  the  3D  wave 
speed  model  allows  prediction  of  complexity  in  body 
wave  arrivals  and  more  pronounced  surface  wave 
arrivals  and  dispersion.  The  CMT  inversion 
procedure  automatically  identifies  windows  to  use  as 
constraints  (Maggi  et  al,  2009),  and  use  of  the  3D 
wave  speed  model  results  in  more  time  windows  and 
subsequently  a  larger  portion  of  the  hill  waveform 
included  in  the  CMT  inversion.  Windows  must 
exceed  stringent  minimums  over  all  metrics  for 
inclusion  in  the  inversion. 

Performance  of  the  CMT  inversion  is  highlighted  in 
Figure  3.  The  assess  the  performance  of  the  CMT 
inversion,  the  p-axis  was  used  as  it  reduces  a 
complex  system  with  up  to  9,  Mij,  parameters  down 
to  2,  strike  and  dip.  This  assumes  the  events  are 
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primarily  double  couple  with  minimal  isotropic  or 
compensated  linear  vector  dipole  components. 

Results  from  a  bootstrap  procedure  on  a  single  event 
are  shown  on  in  the  left  panel  of  Figure  3.  Focal 
spheres  for  ID  and  3D  wave  speed  models  are  shown 
for  reference.  Similar  to  what  is  demonstrated  in 
Figure  3,  for  each  event  and  each  wave  speed  model, 
a  bootstrap  was  performed  to  determine  the  standard 
error  of  the  strike  and  dip  of  the  p-axis.  Standard 
errors  for  the  entire  data  set  for  the  strike  and  dip  are 
shown  in  Figure  3,  center  and  right  panels.  Over  the 
full  data  set,  the  standard  error  of  the  solution 
decreases  when  using  the  3D  model  over  the  ID 
model  at  25-125  seconds  period.  The  reduction  in 
solution  error  is  also  apparent  at  periods  of  60-125 
seconds. 


Figure  2.  Source  database  derived  from  the  3D  wave  Adioint  '“versions  and  Model  Updates 

speed  model  at  25  second  and  longer  t  1  1 

Using  the  differences  between  the  synthetic  and 

observed  seismograms,  an  event  kernel  is  constructed 

for  each  source  using  the  adjoint  method  (Tromp  et  al.  2005).  The  sum  of  all  these  event  kernels  guides  model 
updates  for  the  current  iteration.  Automatic  measurements  are  made  between  the  data  and  synthetics  using  the 
Flexwin  tool  (Maggi  et  al.,  2009),  creating  data  windows  and  associated  metrics  that  are  then  used  in  the  creation  of 
adjoint  sources.  Windows,  and  their  adjoint  sources,  are  not  constrained  to  any  specific  seismic  phase  or  time 
window  and  only  use  well  matched  data- synthetic  pairs.  Adjoint  sources  along  with  the  reconstructed,  full  synthetic 
wavefield  are  propagated  in  reverse  time  through  the  current  iteration  wave  speed  model  to  identify  locations  in  the 
model  requiring  improvement.  Interactions  between  the  adjoint  sources  and  the  time-reversed  wavefield,  integrated 
over  time,  generate  kernels  specific  to  each  measurement,  e.g.  a  simple  difference  between  the  data  and  synthetics. 
Such  kernels  can  be  created  using  individual,  isolated  measurements  or  for  each  event  with  a  large  number  of  back- 
propagated  adjoint  sources.  Event  kernels  are  then  summed  volumetrically  to  produce  Frechet  derivates  that  are  used 
to  update  the  wave  speed  model. 


In  a  first  attempt,  our  updates  to  the  wave  speed  model  used  a  steepest  descent  method  and  an  isotropic  mantle  wave 
speed  along  with  tight  controls  on  the  magnitude  of  the  update.  The  maximum  change  in  the  model  update  was 
limited  to  3%  based  on  comparisons  between  data  and  synthetics  from  a  set  of  test  events  not  used  in  the  inversion. 
As  a  result  of  the  parameterization  the  shear  wave  speed  model  oscillated  between  each  iteration  step.  The  inversion 
was  attempting  to  fit  the  Rayleigh  and  Love  waves  simultaneously.  Due  to  the  long  period  nature  of  these  surface 
waves  they  are  primarily  sensitive  to  mantle  wave  speeds.  An  isotropic  parameterization  of  the  mantle  is  unable  to 
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Figure  3.  Performance  of  the  CMT  inversion  when  using  the  ID  (white)  and  3D  (gray)  wave  speed  models. 
Left  panel  shows  a  histogram  of  p-axis  strikes  from  a  bootstrap  procedure  for  a  single  event. 
Focal  spheres  for  ID  and  3D  with  all  p-axes  from  the  bootstrap  shown  for  reference.  Center 
(right)  panel  shows  the  standard  error  of  the  strike  (dip)  of  the  p-axis  for  the  entire  data  set.  Use 
of  the  3D  wave  speed  model  reduces  the  standard  error  over  the  entire  data  set. 
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Figure  4.  Comparison  of  travel  time  measurements  between  the  initial  model,  gray,  and  the  current 

iteration  model,  colors,  made  using  cross  correlation.  Measurements  are  separated  by  wave  type 
including  Rayleigh  (Rayl)  and  Love  waves,  and  body  waves  on  the  vertical-radial  components 


satisfy  the  known  transversely  isotropic  mantle  and  as  such  the  subsequent  model  iterations  oscillated  between  two 
local  minima.  From  these  initial  tests,  addition  of  an  anisotropy  parameterization,  as  transverse  isotropy,  of  the 
mantle  fit  the  Rayleigh  and  Love  wave  simultaneously  and  allowed  further  incorporation  of  shorter  period  arrivals. 
With  this  new  parameterization,  finer  scale  structures  began  to  appear  and  large-scale  oscillations  were  eliminated 
with  each  iteration. 

Figure  4  shows  the  travel  time  performance  of  the  current  iteration  model  against  the  starting  wave  speed  model.  In 
general,  travel  time  variations  are  reduced  at  long  periods  for  the  Rayleigh  and  Love  waves.  The  total  misfit  between 
the  data  and  synthetics  as  measured  by  metrics  from  Flexwin,  time  shift,  amplitude  ratio,  and  cross-correlation 
value,  all  show  progressive  improvement  to  the  synthetics  relative  to  the  data  for  each  iteration.  Amplitude  ratios 
show  less  of  an  improvement  when  compared  to  travel  time.  The  small  improvement  in  amplitude  is  due  to  a 
combination  of  the  lack  of  even  shorter  periods  in  the  inversion  and  incorrect  attenuation  structure.  Currently  a  ID 
attenuation  structure  is  used  throughout  the  model  domain,  but  addition  and  improvement  to  a  3D  attenuation  model 
are  realizable  by  using  the  amplitude  misfit  for  all  arrivals  and  the  adjoint,  amplitude  kernels  defined  in  Tromp  et  al. 
(2005). 

Progression  of  the  wave  speed  model  for  a  sample  set  of  iterations  is  shown  in  Figure  5  as  absolute  shear  wave 
speeds,  Vsv  and  Vsh.  The  oscillations  previously  identified  are  not  apparent  in  this  set  of  iterations  due  to  the 
incorporation  of  transverse  isotropy,  Vsv  and  Vsh.  As  the  number  of  iterations  increase,  improvements  are  applied 
to  the  model  and  shorter  period  date  are  included  in  the  inversion.  With  these  improvement  and  shorter  period  data 
higher  resolution  structures  begin  to  appear  in  the  wave  speed  model.  Two  examples  of  new  higher  resolution 
structures  are  in  the  Red  Sea  in  Vsv  and  in  the  South  Caspian  Sea  in  Vsh.  It  is  interesting  to  note  the  wave  speed 
model  prefers  reductions  in  wave  speed  in  Vsv  and  increases  in  wave  speed  as  Vsh;  a  feature  that  would  not  be 
possible  without  transverse  isotropy  in  the  mantle  and  the  fitting  of  both  Love  and  Rayleigh  waves. 

A  set  of  comparisons,  waveforms  and  cross-sections,  between  the  current  iteration  and  the  starting  model  (S2.9EA) 
are  contained  in  Figure  6.  Waveform  fits  from  a  path  from  southern  Zagros  to  station  KBK  demonstrates  that 
updates  to  the  model  are  improving  the  fit  between  data  and  synthetics.  This  event  is  not  included  in  the  adjoint 
tomography,  but  is  used  to  provide  a  consistency  check  on  the  models  improvement  and  progression.  After  about  a 
third  of  the  total  iterations  a  lower  wave  speed  Vsv  lithosphere  is  required  in  the  middle  of  this  path,  to  the 
southwest  of  Tibet’s  western  syntaxis.  This  lower  wave  speed  region  appears  in  the  cross  sections,  but  also  apparent 
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in  the  synthetics.  The  surface  wave  in  the  starting  model  arrives  earlier  than  the  data,  but  arrives  “on-time”  using  the 
current  iteration  model. 

The  initial  model  in  central  Asia  at  1 00  km  depth  shows  faster  wave  speeds  caused  primarily  by  thick,  cold 
lithosphere.  Further  to  the  south,  into  the  Zagros  Mountains  where  a  large  amount  of  continental  deformation  and 
crustal  thickening  is  occurring,  the  wave  speeds  are  substantially  slower.  This  Asian  north-to-south,  fast-to-slow 
wave  speed  pattern  in  the  initial  model  is  mirrored  in  the  current  iteration,  Figures  5  and  6.  Similar  behavior  was 
identified  in  other  finite-frequency  tomography  experiments  (Montelli  et  al.,  2004).  Wave  speed  perturbations  in  the 
initial  model  are  amplified  when  using  sensitivity  kernels  based  on  a  finite-frequency  approach  rather  than  an 
approximation,  either  classical  ray  theory  or  surface  wave  sensitivities  kernels.  Large  continuous  and  anomalous 
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Figure  5.  Comparison  of  the  Middle  East  wave  speed  model  at  different  iterations.  The  left  column  is 
Vsv  and  the  right  column  is  Vsh.  The  top  row  is  the  initial  wave  speed  model,  the  second  row 
is  the  3rd  iteration  and  the  bottom  row  is  the  6th  iteration.  Model  depth  is  100  km. 
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regions  in  the  initial  model,  Figure  5,  will  become  stronger  and  begin  to  separate  into  individual  anomalies  as  the 
more  iterations  are  completed. 

Updates  to  the  model  also  include  density  and  compressional  components.  The  initial  3D  compressional  wave  speed 
model  was  scaled  in  the  initial  model  from  the  3D  shear  wave  speed  by  a  constant  factor.  Shear  wave  speeds  in 
S2.9EA  were  derived  from  a  large  number  of  Love,  Rayleigh,  and  shear  body.  A  comparison  between  the 
compressional  adjoint  results  and  the  initial  model  can  be  difficult  as  the  ratio  between  shear  and  compressional 
wave  speeds  can  vary  due  to  composition.  All  adjoint  iterations  for  compressional  wave  speeds  show  a  much  slower 
upper  mantle  beneath  much  of  the  Middle  East  than  indicated  by  S2.9EA.  Updates  to  the  model  at  long  wavelengths 
for  these  initial  iterations  are  predominately  controlled  by  Rayleigh  waves  (shear  and  compressional  surface  wave) 
as  the  initial  S2.9EA  model  fit  the  Love  waves  (shear  only  surface  waves)  better,  Figures  4  and  5. 


Figure  6.  Demonstration  of  the  current  iteration  model  through  a  single  station  -  event  pair.  The  station, 
KBK,  and  event,  southern  Zagros,  are  not  included  in  the  adjoint  tomography.  Panel  a)  includes 
a  location  map  and  a  cross-section  with  relative  wave  speeds  between  the  initial  model,  mO,  and 
the  current  model,  m.  Panel  a)  also  includes  the  same  cross-section  but  as  absolute  wave  speed, 
Vsv.  Panel  b)  shows  the  improved  fit  of  the  waveforms  between  the  starting,  left,  and  current 
models,  right.  Data  is  black  and  synthetics  are  plotted  in  color. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Adjoint  iterations  to  update  the  wave  speed  model  of  the  Middle  East  have  shown  dramatic  improvements  to  the 
model  based  on  data  /  synthetic  comparisons  and  expected  behavior  when  introducing  a  finite-frequency 
tomographic  approach.  Before  iterations  began,  a  seismic  waveform  database  with  -200  reinterpreted  sources  in  the 
Middle  East  was  compiled  to  avoid  the  mapping  of  source  errors  into  wave  speed  structure.  All  sources  in  the 
database  agree  with  previously  published  solutions  and  between  the  methodologies  used:  CAP  (Zhao  and 
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Helmberger,  1994),  teleseismic  (Kikuchi  and  Kanamori,  1982)  and  a  CMT  inversion  (Liu  et  al.,  2004).  Initial 
comparisons  between  data  and  synthetics  computed  from  the  initial  3D  model  S2.9EA  (Kustowski  et  al.,  2008)  at  25 
seconds  and  longer  show  reasonable  agreement  and  a  general  improvement  over  the  ID  PREM  model  (Dziewonski 
and  Anderson,  1981).  To  further  improve  the  fit  between  data  and  synthetics,  the  adjoint  inversion  methodology  was 
implemented  to  iteratively  update  the  3D  wave  speed  model  (Tromp  et  al.,  2005).  Following  an  initial  set  of 
iterations,  the  wave  speed  model  was  parameterized  to  incorporate  transverse  isotropy  in  the  mantle  to  accommodate 
the  simultaneous  fitting  of  Rayleigh  and  Love  waves.  This  parameterization  helped  to  stabilize  the  inversion  and 
allowed  further  inclusion  of  shorter  periods  into  the  adjoint  tomography.  As  the  model  progresses  through  a  series  of 
iterations,  regions  of  wave  speed  anomalies  have  increased  in  strength  from  the  initial  to  the  current  iteration  as 
physics  of  wave  propagation  were  incorporated  into  the  inversion;  see,  e.g.,  the  Red  Sea  and  the  southern  Caspian 
Sea,  Figure  5.  Compressional  wave  speed  estimates  also  improved  from  increased  use  of  Rayleigh  waves.  Our 
recommendations  for  a  robustly  constrained  3D  wave  speed  model  in  the  Middle  East  include  using  quality 
controlled  seismic  waveform  data  from  well-constrained  sources  and  use  of  a  physics  based  tomographic  technique 
to  improve  upon  existing  3D  models.  Further  recommendations  to  improve  the  wave  speed  model  include  adding: 

•  anisotropic  effects  where  necessary,  e.g.,  the  mantle,  and 

•  attenuation  to  improve  fits  to  amplitudes,  and 

•  shorter  period  signals  to  further  improve  the  waveform  predictive  capabilities  of  the  wave  speed  model  of 
the  Middle  East. 
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